####################################################
### Estimate Demand Model for Covered California ###
####################################################
		
##### Specify variables and Estimation Procedure

	# ssubmit --partition=sscc --cores=3 --mem=200g "julia estimate.demand.inertia.jl"
	# ssubmit --partition=sscc --cores=35 --mem=150g "julia18 estimate.demand.inertia.jl"
	cluster_run = false;

	eq_robust_flag = false;
	network_flag = false;
	use_newton = true; # Nelder-Mead if false	
	control_function = true;
	random_parameters = true;
	sequence_flag = false;
	inattention = false;
	nested = false;
	
	add_premium_interactions = true;
	add_previous_demographics_interactions = true;
	add_previous_firm_interactions = true;
	add_av_interactions = true;
	add_av_interactions_noinc = false;
	add_intercepts = true;
	add_rating_areas = true;
	add_silver = false;
	add_insurer_market_fe = true;
	add_complex_interactions = true;
	add_random_cov_terms = false;
	add_av_as_lognormal = false;
		
	shrink_sample = false;
	smooth_flag = false;
	do_estimation_flag = false;
	tau = 1;	
		
	year_to_run = 2018; # supply model uses data from 2014-2018
	
	use_eta = false;
	sequence_exp_flag = false;
	
	keep_years_in_panel = 0;
	keep_two_choices = false;
	keep_three_choices = false;
	sample_size = 0.20;
	save_covariate_flag = false;
	ns = 50;
	
	rich_flag = false;
	poor_flag = false;
	old_flag = false;
	young_flag = false;
	lapse_flag = false;
	no_lapse_flag = false;
	churn_flag = false;
	no_churn_flag = false;
	
	if network_flag
		fill_missing = false;
		keep_missing = false;
		network_interaction_flag = false;
	end
	
	if cluster_run
		codedir = "/project/inertia_project/Code";
		datadir = "/project/inertia_project/Data";
	else
		codedir = "C:/Users/eas24f/OneDrive - Florida State University/Research/Inertia RESTAT/Code"; 
		datadir = "C:/Users/eas24f/OneDrive - Florida State University/Research/Inertia RESTAT/Data"; 
	end
	
##### Specify starting point and solution method
	
	cd(datadir)
	using CSV
	using DataFrames
	starting_point = CSV.read("starting_point.csv",DataFrame,header = true); # Read from file - modify file to change starting point	
	beta_initial = starting_point[!,:Coef];		
	if smooth_flag
		true_solution = starting_point[!,:Truth];		
	end
	
##### Load Support Functions and Data and Setup Covariate Matrix
	
	cd(codedir)
	include("demand.functions.inertia.jl");
	include("elasticity.functions.CA.jl");
	include("setup.demand.data.jl");
	X = form_X_matrix(all_covariates);	
	if inattention
		Xsearch = form_Xsearch_matrix(search_covariates);
	end
	
##### Estimation
	
	# Run using Nelder Mead or Newton's
	if nested & do_estimation_flag
		if use_newton
			min_beta = newton_nested_logit(beta_initial,X,choices,1e-8,control_function);
		else
			obj_nested_logit(betap) = objective_nested_logit(betap,X,choices,control_function,random_parameters);
			nested_logit = optimize(obj_nested_logit,convert(Array,beta_initial),iterations=50000,show_trace=true) ;
			min_beta = Optim.minimizer(nested_logit)
		end
	elseif do_estimation_flag
		if use_newton
			min_beta = newton_logit(beta_initial,X,choices,1e-8,control_function);
		else
			obj_logit(betap) = objective_logit(betap,X,choices,control_function,random_parameters);
			logit_results = optimize(obj_logit,convert(Array,beta_initial),iterations=50000,show_trace=true) ;
			min_beta = Optim.minimizer(logit_results)
		end
	else
		min_beta = copy(beta_initial);
	end
	
	probability_matrix = get_choice_probabilities(min_beta);
	if inattention
		prob_matrix = copy(probability_matrix);
		probabilities = vec(sum(dims=2,probability_matrix));
	elseif random_parameters & sequence_flag
		prob_matrix = probability_matrix[1:length(choices),:];
		sequence_prob_matrix = probability_matrix[(length(choices)+1):size(probability_matrix)[1],:];
		probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
		sequence_probs = vec(sum(dims=2,sequence_prob_matrix)) ./ ns;
	elseif random_parameters
		prob_matrix = copy(probability_matrix);
		probabilities = vec(sum(dims=2,prob_matrix)) ./ ns;
	else
		probabilities = probability_matrix;
	end
	
	# Simulation exercise to check random parameters estimation is working
		
	check_random_parameters_flag = false;
	if check_random_parameters_flag
		logit = check_random_parameters(beta_initial,true_solution);
	end
	
##### Output Processing

	save_output = true; # save estimates
	if save_output
		param_output = process_parameters(min_beta,probabilities);
		CSV.write(filename,param_output);
	end
	
	save_choice_output_flag = true;
	if save_choice_output_flag
		choice_output = DataFrame();
		choice_output[!,:prob] = probabilities;
		if control_function
			choice_output[!,:instrumented_prob] = calculate_instrumented_probs(min_beta);
		end
		CSV.write(choice_output_filename,choice_output);
	end
	
	compute_elasticities_flag = true;
	if compute_elasticities_flag
		elasticity_output = compute_elasticities(min_beta,probabilities);
		CSV.write(elasticities_filename,elasticity_output);
	end

	compute_switching_costs_flag = true;
	if (compute_switching_costs_flag) & (!inattention)
		inertia_output = compute_switching_costs(min_beta,probabilities);
		CSV.write(inertia_output_filename,inertia_output);
	end
	
	compute_average_params_flag = true;
	if compute_average_params_flag
		avg_param_output = compute_average_params(min_beta,probabilities)
		CSV.write(avg_param_output_filname,avg_param_output);
	end
	
	